{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Area under 1-D and 2-D curves, various methods\n",
"\n",
"**Randall Romero Aguilar, PhD**\n",
"\n",
"This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.\n",
"\n",
"Original (Matlab) CompEcon file: **demqua03.m**\n",
"\n",
"Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
"\n",
" !pip install compecon --upgrade\n",
"\n",
"Last updated: 2022-Oct-23\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## About\n",
"\n",
"Uni- and bi-vaiariate integration using Newton-Cotes, Gaussian, Monte Carlo, and quasi-Monte Carlo quadrature methods."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Initial tasks"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from compecon import qnwtrap, qnwsimp, qnwlege\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"quadmethods = [qnwtrap, qnwsimp, qnwlege]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make support function"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"a, b = -1, 1\n",
"nlist = [5, 11, 21, 31]\n",
"N = len(nlist)\n",
"\n",
"def quad(func, qnw, n):\n",
" xi, wi = qnw(n,a,b)\n",
" return np.dot(func(xi),wi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Evaluating\n",
"$\\int_{-1}^1e^{-x}dx$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.exp(-x)\n",
"\n",
"f_quad = np.array([[quad(f, qnw, ni) for qnw in quadmethods] for ni in nlist])\n",
"f_true = np.exp(1) - 1/np.exp(1)\n",
"f_error = np.log10(np.abs(f_quad/f_true - 1))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Evaluating\n",
"$\\int_{-1}^1\\sqrt{|x|}dx$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def g(x):\n",
" return np.sqrt(np.abs(x))\n",
"\n",
"g_quad = np.array([[quad(g, qnw, ni) for qnw in quadmethods] for ni in nlist])\n",
"g_true = 4/3\n",
"g_error = np.log10(np.abs(g_quad/g_true - 1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Make table with results"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" Trapezoid rule | \n",
" Simpson's rule | \n",
" Gauss-Legendre | \n",
"
\n",
" \n",
" \n",
" \n",
" | $\\int_{-1}^1e^{-x}dx$ | \n",
" 5 | \n",
" -1.683044 | \n",
" -3.472173 | \n",
" -9.454795 | \n",
"
\n",
" \n",
" | 11 | \n",
" -2.477411 | \n",
" -5.053217 | \n",
" -14.273349 | \n",
"
\n",
" \n",
" | 21 | \n",
" -3.079254 | \n",
" -6.255789 | \n",
" -14.675836 | \n",
"
\n",
" \n",
" | 31 | \n",
" -3.431396 | \n",
" -6.959867 | \n",
" -15.653560 | \n",
"
\n",
" \n",
" | $\\int_{-1}^1\\sqrt{|x|}dx$ | \n",
" 5 | \n",
" -1.023788 | \n",
" -1.367611 | \n",
" -0.870112 | \n",
"
\n",
" \n",
" | 11 | \n",
" -1.595301 | \n",
" -1.347900 | \n",
" -1.351241 | \n",
"
\n",
" \n",
" | 21 | \n",
" -2.034517 | \n",
" -2.414470 | \n",
" -1.758970 | \n",
"
\n",
" \n",
" | 31 | \n",
" -2.293296 | \n",
" -2.063539 | \n",
" -2.007803 | \n",
"
\n",
" \n",
"
\n",
"